%	Compute the comparative statics of risk measures with respect to
%	changes of beta=delta

clear all; close all;
tic
warning('off')

cd ..\
cd ..\
cd Matlab_functions

global wgrid betta delta gamma pi1 pi2 theta e1 e2 ey1 ey2 n v1aut v2aut wfb vfb nw w_st wmin wmax

%----------------------------------------------
% Parameter Setting
%----------------------------------------------

gamma = 1;                              % Risk Aversion
pi1 = 0.5;                              % prob state 1
pi2 = 1-pi1;                            % prob state 2
theta = 0;                              % aggregate shock  
e1 = 1+theta;                           % endowment state 1
e2 = 1-(pi1/(pi2))*theta;               % endowment state 2
kappa = 0.6;
sigma= 0.1;                 % measure of mean of endowment
ey1 = (kappa-sigma*pi2/pi1)*e1;         % endowment state 1 young
eo1 = (1-kappa+sigma*pi2/pi1)*e1;       % endowment state 1 old
ey2 = (kappa+sigma)*e2;                 % endowment state 2 young
eo2 = (1-kappa-sigma)*e2;               % endowment state 2 old                
n=0;                                    % fertility rate
delta_vector=[0.65,0.70,0.75,0.80,0.9,0.95,exp(-1/75)];   

for hh=1:length(delta_vector)
  
delta=delta_vector(hh);
betta=delta;

%----------------------------------------------
% Checks
%----------------------------------------------

lamb1 = (eo1/ey1)^gamma;
lamb2 = (eo2/ey2)^gamma;
S1 = lamb1-betta;          % this must be positive, state 1 bad state for the young
S2 = lamb2-betta ;         % this must be negative, state 2 good state for the young
Sfb = lamb1-betta/delta;   % if negative, then FB implies positive transfers in both states
No_Samuelson = (1-kappa)-betta*kappa;   % if positive, then dynamic efficiency in deterministic case

%----------------------------------------------
% Find max and min of promised values
%----------------------------------------------

[v1aut,v2aut,waut] = Autarky(gamma,betta,pi1,pi2,ey1,ey2,eo1,eo2);
w_st = Maxstat(gamma,betta,pi1,pi2,ey1,ey2,eo1,eo2,v1aut,v2aut);

%----------------------------------------------
% Find the maximum interior convergence
%----------------------------------------------

[vfb,wfb,cy1fb,cy2fb,w_conv]=Intconv(gamma,betta,delta,n,pi1,pi2,v1aut);

%----------------------------------------------
% grid for promised values
%----------------------------------------------

wmin = waut;
wmax = w_st;
nw = 300; % search method for optimal grid number
wgrid = cheby_points(wmin,wmax,nw)';

%----------------------------------------------
% Iterating the Policy Function
%----------------------------------------------

distmia=100;
olddist=1;
iter=0;
tollv=10^(-4);
tic
fprintf('For this tolerance, I predict that about %g iterations will be needed \n',round(fsolve(@(x)tollv.^(1/x)-delta,3)));

cy1_old = cy1fb*ones(nw,1);
cy2_old = cy2fb*ones(nw,1);
w1_old = wfb*ones(nw,1);
w2_old = wfb*ones(nw,1);

% Policy iteration 
P_old=[cy1_old,cy2_old,w1_old,w2_old];
 
 while distmia > tollv
    iter=iter+1;
 
    [V,cy1,cy2,w1,w2,w_thr,w_approx] = Value_Iteration(cy1_old,cy2_old,w1_old,w2_old); 
    w_upper=fsolve(@(x)(x-interp1(wgrid,w2,x,'spline')),waut);

    P=[cy1,cy2,w1,w2]; 
    
    distmia=max(max(abs(P-P_old)));
    distmia=distmia+abs(cy1fb-interp1(wgrid,cy1,w_thr,'spline'));
    
    fprintf('Iteration %g completed, Distance is: %g   Improved by %9.1f %% \n',iter,distmia,(((olddist-distmia)/olddist))*100);
     
    olddist=distmia;
    P_old=P;
    cy1_old=cy1;
    cy2_old=cy2;
    w1_old=w1;
    w2_old=w2; 
 end
 
kk=200;
w_thr=w1(1);
wmax = fsolve(@(x)interp1(wgrid,w2,x,'spline')-x,w1(1));
gridErg = [w_thr:(wmax-w_thr)/kk:wmax];
nw = length(gridErg);

w1 = (interp1(wgrid,w1,gridErg,'spline'));
w2 = (interp1(wgrid,w2,gridErg,'spline'));
cy1 = (interp1(wgrid,cy1,gridErg,'spline'));
cy2 = (interp1(wgrid,cy2,gridErg,'spline')); 

% next period consumption
cy11=interp1(gridErg,cy1,w1,'spline');
cy12=interp1(gridErg,cy2,w1,'spline');
cy21=interp1(gridErg,cy1,w2,'spline');
cy22=interp1(gridErg,cy2,w2,'spline');

% Kernel prices (Stochastic discount factors)
m11 = betta.*(cy1./(e1-interp1(gridErg,cy1,w1,'spline'))).^gamma;
m12 = betta.*(cy1./(e2-interp1(gridErg,cy2,w1,'spline'))).^gamma;
m21 = betta.*(cy2./(e1-interp1(gridErg,cy1,w2,'spline'))).^gamma;
m22 = betta.*(cy2./(e2-interp1(gridErg,cy2,w2,'spline'))).^gamma;

V0 = (interp1(wgrid,V,gridErg,'spline'));

% consumption equivalent welfare
CE=(exp(((delta.*(1-delta))./(betta+delta)).*(vfb-V0))-1).*100;

% Insurance coefficient
IC1=1-(log(cy11)-log(cy12))./(log(ey1)-log(ey2));
IC2=1-(log(cy21)-log(cy22))./(log(ey1)-log(ey2));

% transition matrix
edge= gridErg(1:nw);
edge(nw+1)=gridErg(nw);
gridk = discretize(gridErg,edge);
w1kk = discretize(w1,edge);
w2kk = discretize(w2,edge);

w1kk(1)=1;

for i=2:nw
if isnan(w2kk(i))
   w2kk(i)=w2kk(i-1); 
    
end

if isnan(w1kk(i))
   w1kk(i)=w1kk(i-1); 
    
end

end

[T,M,Q]=Matrices(m11,m12,m21,m22,w1kk,w2kk,gridk,nw,pi1);
T=T';

% stationary distribution
probw = (1/(2*nw))*ones(2*nw,1);
   test=1;
   while test > 12^(-12);
      probwx = T*probw;
      test = max(abs(probwx-probw));
      probw = probwx;
   end;

probw=probw/sum(probw); 

for i=1:nw
  probw1(i)=probw(i);
end

sum(probw1)

for i=nw+1:2*nw
   probw2(i-nw)=probw(i); 
end
sum(probw2)

probw = probw1+probw2;

%Iterative method to find Ross measure
x0 = ones(length(Q),1);

distmia=100;
iter=0;
tollv=10^(-10);

while distmia>=tollv
iter = iter+1;

x1 = Q*x0;
yy0 = x0./max(x0);
yy1 = x1./max(x1);

distmia = max(abs(((yy1)-(yy0))./((yy0))));

x0 = x1;
end

max_eig = yy1;

% Ross measure of risk
Ross_betadelta(hh) = -log(min(max_eig)/max(max_eig))

% Consumption equivalent welfare change
CE_betadelta(hh)=probw*CE'

% Insurance coefficient
IC_betadelta(hh)=pi1*probw*IC1'+(1-pi1)*probw*IC2'

end
 
cd ..\
cd Stored_File

save('Comparative_betadelta.mat','delta_vector','Ross_betadelta','IC_betadelta', 'CE_betadelta')
